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ABSTRACT 



An investigation was ma 




uckling of rectangular plates 



with free, simply supported, or clamped edges and loaded in uniform 
compression or pure shear. The finite element method with a rectangu- 
lar, sixteen-degree-of-freedom element was used. Results show that, 
for the same order of error, this element allows much coarser plate 
divisions than a widely used twelve-degree-of -freedom element. An 
order of magnitude reduction in the required order of the assembled 
stiffness and stability coefficient matrices results. 
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LIST OF SYMBOLS 



a plate length in x direction; polynomial coefficient 

a vector of displacement polynomial coefficients 

A transformation matrix 

b plate length in y direction 

B stability coefficient matrix 

D plate flexural rigidity, Eh^/(12(1 - ^ ^)) 

E modulus of elasticity 

i vector of external nodal forces and moments 

G matrix of functions of x and y, relating vectors a 

and z 

h plate thickness 

H matrix of functions of x and y, relating vectors a 

and r 

I. identity matrix 

k buckling coefficient 

K stiffness matrix 

K* modified stiffness matrix from defining relation K - DK* 
m vector of functions of bending moments 

M bending moment per unit length 

N* load matrix 

N matrix of constants giving relative magnitudes of types 

of loads 

p plate element length in x direction 

P matrix of functions of Poisson D s ratio 

q plate element width in y direction 

£ vector of first derivatives of transverse displacement 

U strain energy of bending 

V total potential energy 
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1. INTRODUCTION 



BACKGROUND 

The Finite Element Method for Stability Problems 

The finite element method provides a numerical means of analyzing 
boundary value field problems * When applied to structural mechanics, 
the method provides a means of deriving a stiffness matrix for an as- 
sembly of structural elements. For stability problems, a stability co- 
efficient matrix is also required. This matrix is used in conjunction 
with the stiffness matrix. 

Development of the Method 

Turner, Clough, Martin, and Topp introduced the finite element 
1 

method in 1956. In 1963, Melosh presented a basis for derivation of 

the stiffness matrix for a thin plate element, using three degrees of 
2 

freedom per node. Also in 1963, Gallagher and Padlog first used a 

3 

stability coefficient matrix, for column buckling. In 1966, Kapur and 

Hartz gave the results of the use of the stability coefficient matrix 

4 

approach for plate buckling problems. They used the Melosh plate ele- 
ment to find the stiffness matrix. While Kapur and Hartz were making 

their study of plate buckling, in 1965 Bogner, Fox, and Schmit presented 

5 

a new stiffness matrix for a plate element. Four degrees of freedom 
per node were used for the new element. 

Significance of Degrees of Freedom 

For the Melosh method of deriving the stiffness matrix of a plate 
element, the three degrees of freedom per node are lateral deflection 
and slope in each of two perpendicular, in-plane directions (w; w,^; 
w,^). This assures compatibility of lateral deflection of adjacent edges 
of elements. However, complete continuity of slope is not achieved. The 
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method of Bogner, Fox, and Schmit includes an additional degree of 



freedom (w, ) which does provide continuity of slope. For a rectangu- 

xy 

lar plate element, an element stiffness matrix of order sixteen is re- 
quired to represent four degrees of freedom at each of four nodes. For 
three degrees of freedom per node, the order is twelve. Hence, the six- 
teen-degree-of-freedom element stiffness matrix is larger but should 
produce more accurate results because of better continuity of slope.* 



THE PROBLEM INVESTIGATED 

Ob jective 

The objective of the investigation was to determine the consequences 
of the use of the newer, sixteen-degree-of -freedom element stiffness 
matrix in plate buckling problems. Secondarily, some previously unsolved 
cases of plate buckling were investigated to show the versatility of the 
finite element method. 

Limitations 

The investigation was limited to finding buckling coefficients and 
buckling mode vectors for thin, homogeneous, isotropic, rectangular plates 
with free, simply supported, or clamped edge conditions, and with uniform 
compression or pure shear loads. Plates were subdivided into identical 
rectangular elements. 

JUSTIFICATION OF THE INVESTIGATION 
Advantages of the Finite Element Method 



*Other rectangular plate elements with twelve or sixteen degrees of 
freedom may be formed by using different displacement functions. Through- 
out, plate elements with twelve or sixteen degrees of freedom will mean 
the elements discussed above. 
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Plate buckling problems of even moderate complexity require numeri- 
cal procedures o Of the numerical procedures applicable, the finite ele- 
ment method has some unique advantages, as follows: 

lo The capability of analyzing plates of arbitrary shape, such 
as plates with holes, 

2, The capability of analyzing plates made of materials with 
orthotropic properties, 

3o The capability of introducing edge conditions conveniently, 
especially the free edge, 

4. The capability of synthesizing complex structures from plates 
and other members, such as beams. 

Advantages of Using Sixteen-Degree-of-Freedom Element 

Assembled stiffness and stability coefficient matrices may be very 
large, especially in structural problems which require synthesis of 
several types of members. Finer divisions of a plate require larger 
assembled matrices. Any procedure which would reduce the order of 
these matrices, without increasing the error of the result, would, of 
course, save computer time and space. 

It was anticipated that use of the newer element would allow coarser 
plate divisions with no increase of error. If so, a significant reduc- 
tion in the size of the assembled matrices might be possible in spite of 
the higher order of the element matrices. Hence, the investigation was 
undertaken to see whether the newer element produces a significant im- 
provement , 



PROCEDURE 

The equations for bending energy and work of in-plane forces were 
written in matrix form, making use of the finite element method with 
the sixteen-degree-of-f reedom element, A computer program, with matrix 
iteration, was written to generate and solve the matrix equation re- 
sulting from application of conditions of minimum potential energy. 
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Several previously solved cases of plate buckling were run. Results 
were compared with known values and with values obtained using the twelve- 
degree-of-f reedom plate element. A few previously unsolved cases were 
investigated. 
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2. FORMULATION 



INTRODUCTION 

In this section, the formulation of the problem is introduced by 
a discussion of the energy method* The plate element, displacement 
vector, displacement polynomial, and transformation matrix are explained 
as preliminaries. With this background established, the matrix equa- 
tions for bending energy and work of in~plane forces are derived for a 
single plate element. Extension of the method to multiple-element 
plates and the method of handling boundary conditions are discussed next. 
Then, the method of finding a matrix equation for minimum potential energy 
is shown* Finally, the solution of the resulting characteristic value 
problem is discussed. 



THE ENERGY METHOD 

If a plate loaded by external, in-plane forces is bent a small 
amount, the external forces move because the plate edges move. Thus, 
the forces do work. At the same time, to bend the plate, the strain 
energy of bending must be supplied. If the work of external forces 
exceeds the bending energy, the plate buckles. These two energy terms 
are equal in magnitude at the critical load. Or, the total potential 
energy is minimum at the critical load. 

PRELIMINARY DEFINITIONS 
Plate Element and Displacement Vector 
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(w) 




The plate element used is shown above. The nodal order used is 
indicated at the corners of the element. The generalized displacement 
vector used for a single node is* 

tfl = C w w, y w,„ w, Ky :i T (1) 

The displacement vector used for a plate element is 



An underline indicates a matrix. Capital letters are used for 
square or rectangular matrices. Small letters are used for column 
vectors. A superscript T is used to indicate a transpose. 
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rfo = C rfi(0,0) i c£,(p,0> ; tf 3 (0,q) ; tf 4 (p,q^ T 



( 2 ) 



Displacement Polynomial and Transformation Matrix 
Transverse deflection was represented by 

w =Z I a t jx l -y- 1 (3) 

l=i } = i J 

By using this polynomial and derivatives of it, evaluated at the plate 
element nodes, the generalized deflection of the nodes may be represented 
by 



do = Aa (4) 

where A is a 16 x 16 transformation matrix, and a is a vector of the poly- 
nomial coefficients. The polynomial coefficients are then given by 



g = 



(5) 



BENDING ENERGY 

The integral expression for the strain energy of bending is 
U = - |Xf (MxW, x » + M y W;yy - 2. M xy w, X y )dx dy 

This expression may be re-stated in matrix form as 

XJ = Jj m T z. dx dy 



( 6 ) 



(7) 
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where m and z are vectors defined by 

m = L Mx My -ZMxy3 T 



( 8 ) 

(9) 



Z — W»yy W; xy 3^" 

Bending moments are related to curvatures and material properties by 

Mx — — D (Wjxx +■ a) W;yy) 

My = “ D ( W,yy + a) W, xx ) 

Mxy ~ O (1 ~ V ) W, xy 

Here, D is the plate flexural rigidity, and -0 is Poisson's ratio 
Equations (10) may also be written as 

— 

W y xx 



( 10 ) 



Mx 




My 


= - D 


-2M*y 





1 

4 



l 



0 

0 



w. 






(11) 



Defining the matrix of functions of 0 to be P allows equation (11) 
to be given in the more compact form 



m = - QPz 



( 12 ) 



Combining equation (7) and equation (12) gives the new equation for 
bending energy 



U = iDjJ(Pz) T z dxdy 



(13) 



Since the vector z is made up of derivatives of the displacement poly- 
nomial, it may be related to the polynomial coefficients by 



Z = 6a 



(14) 



where G is a 3 x 16 matrix of functions of x and y. Using equation (14) 
to substitute for z_ in equation (13) gives 



TJ = | Dj/(P6a) T 6g dx dy 



(15) 
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Substituting for the vector a using equation (5), extracting terms not 
containing x and y from the integrand, and rearranging gives 

U = | O £ (A' l ) T (JJ G t P G dx dy) A 1 d 0 (16) 

The portion 

K 0 = D(A‘ l ) r (Jj G t P t G dx dy ) A' 1 (17) 

is the stiffness matrix for the plate element. Finally, the bending 
energy of the plate element, in terms of the element stiffness matrix, 
is 

U = | 6j K 0 d 0 (18) 

WORK OF IN-PLANE FORCES 

The integral expression for the additional work done by in-plane 
forces because of lateral deflection is 

W = - 4 f( ( N‘xW,x + Ny W,y +2 Mxy W v x w, y ) dx dy (19) 

The matrix form of this equation is 

W * - |JJ r N* r dx dy (20) 



Here, £ and N* are defined by 



r = 



W,y 

W,x_ 






Ny Nxy 
Nxy Nx 



( 21 ) 



( 22 ) 
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The matrix of load parameters may be expressed as 



= r\N (23) 

where is a scalar factor, and N is a matrix of constants which give 

the relative magnitudes of the types of loads. The vector r, which con- 
sists of derivatives of the displacement polynomial, may be related to 
the polynomial coefficients by 

r= Hd (24) 

H is a 2 x 16 matrix of functions of x and y. Equations (20), (23), 
and (24) lead to 

W = - I j|(H<3) T NBci dxdy (25) 

Continuing in the same manner as was used for bending energy formulation 
gives 

W = - 1 (A' 1 ) T (JjH T NI-j dxdy)A' l tf 0 (26) 

The portion 

B 0 = (A- l ) T (JJ H t N H dx dy ) A -1 (27) 

is the stability coefficient matrix for the plate element.* 

The work of in-plane forces, in terms of the stability coefficient matrix, 
is 

W= - | B 0 <3o (28) 

*The stability coefficient matrix is sometimes referred to as the 
geometric stiffness matrix.® 
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EXTENSION TO MULTIPLE ELEMENTS 



So far, the stiffness matrix and stability coefficient matrix for 
a single plate element have been derived. Poor accuracy results when a 
plate is represented by a single element. Rather, the plate must be 
represented by an assembly of plate elements. Although the plate need 
not be divided into identical elements, identical elements were used 
throughout this investigation. 

The same method of extension to multiple elements for the stiffness 
matrix applies to the stability coefficient matrix. Only a brief dis- 
cussion will be given. A more complete description is contained in 
references 7 and 8. 

Let be a generalized force vector representing externally ap- 
plied forces and moments at nodes of a plate element. Then, and the 
corresponding nodal displacement vector are related by 

fo = Ko < 5 0 ( 29) 

where K q has an order equal to the number of degrees of freedom of the 
plate element. Similarly, when several elements are assembled to make 
up a plate, 

£ = K<5 (30) 

where K is the assembled stiffness matrix, of higher order than K be- 
— — o 

cause there are more nodes. Elements of K are sums of elements taken 
from the plate element stiffness matrices. 

BOUNDARY CONDITIONS 

One of the appealing characteristics of the finite element method is 
the simplicity of handling boundary conditions. For example, no action 
whatever is required for a free edge. For other edge conditions, two ap- 
proaches are available. One method is to eliminate appropriate rows and 

columns of the stiffness and stability matrices, reducing the order of the 
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matrices. For ease of computer programming, a second method was used. 
Zeros were placed in appropriate rows and columns of the matrices. Then 
each zero diagonal element of the stiffness matrix was made unity, be- 
cause the inverse of the matrix was needed later. Table I shows, for 
edge nodes, the nodal displacement components which are zero for the 
prescribed edge condition. It should be noted that the constraints given 
in Table I are the only edge constraints. There is no constraint on in- 
plane deflection for any edge condition; i.e., the plate is "supported” 
or "clamped", but not "held". 

MINIMUM POTENTIAL ENERGY 

Let K and B now represent the assembled stiffness and stability 
matrices, altered by edge conditions. The total potential energy of 
the plate is the sum of bending energy and work of in-plane forces. 

From equations (18) and (28), 

y= U + W = zi T K6 - |r^ T Bcf (31) 

The critical condition for stability is an equilibrium configuration 
with a set of displacements which makes the total energy a minimum. This 
requires that 



K<5 - = O 



(32) 



CHARACTERISTIC VALUE PROBLEM 

With some alteration of form, equation (32) expresses a well-known 

9 

characteristic value problem. 



(DK*)6 = ryB 6 



(33) 
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TABLE I 



BOUNDARY CONDITIONS 



Edge Condition 


w 


w, 

X 


w, 

y 


W, 

xy 


Simply supported edge 
parallel to x axis 


0 


0 






Simply supported edge 
parallel to y axis 


0 




0 




Clamped edge 


0 


0 


0 


0 


Line of symmetry 
parallel to x axis 






0 


0 


Line of symmetry 
parallel to y axis 




0 




0 


Line of anti-symmetry 
parallel to x axis 


0 


0 






Line of anti-symmetry 
parallel to y axis 


0 




0 
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(34) 



(D/^) KM = &6 

A KM = b£ (35)* 

Here, A is an eigenvalue parameter which may have as many values as 
the order of the matrices. Since only the lowest load is desired, 
however, only the largest eigenvalue needs to be found. 

£ 

Most plate buckling coefficients are expressed in the form 

k =(N crit b*)/(ir 2 D) < 36) 

where b is the plate width in the y direction. If the largest value 
of X is known, the buckling coefficient may be given in the same 
form by 

k =(rYb 2 )/('Tr 2 D) = tf/CAir 2 ) (37) 

SOLUTION OF CHARACTERISTIC VALUE PROBLEM 
Solving equation (35) for the largest eigenvalue is the only 
problem remaining. Matrix iteration was the method used, although 
there are other methods. 

The first procedure tried was as follows: 

1. Assume a starting trial vector, 5 . 

2. Calculate an improved trial vector, d', using 

6 ' = K" 1 Bd (38) 

3. Use the Rayleigh quotient to find an estimate for A , using 

A= ((d') T B (39) 

4. Compute k, using 



*This equation is often given in the form (( K* = 0 . 
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k = bV(T^) 



(40) 



5. Use the improved trial vector as the next starting trial 
vector, and iterate until successive values of k agree 
within a prescribed limit. 

The procedure did not give good results for shear loads. For pure shear, 
corresponding to each positive eigenvalue, there is a negative eigen- 
value of equal magnitude. The Rayleigh quotient does not give a good 
estimate for "A . 

A new procedure was used which proved to be satisfactory. A new 
quotient, 

'X 1 =(W) t K6')/(6 t K<5) (4i) 

was used instead of the Rayleigh quotient. Otherwise, the procedure 
was the same as before. 



23 



3. RESULTS 



INTRODUCTION 

General 

In this section, results of three kinds are given. First, compari- 
sons are made with known solutions to show the accuracy of the method. 
Second, comparisons are made with results given by Kapur and Hartz to 
show the effect of the sixteen-degree-of-freedom element. Third, some 
previously unsolved cases of buckling are given to show the versatility 
of the finite element method. 

Information for Reading Tables of Results 

Wherever appearing, the known values of the buckling coefficient 
are from the methods given by Timoshenko and Gere, but are not necessar- 
ily the values stated by them.^ Some values were computed to a greater 
precision. In all cases, the buckling coefficient is given in the form 

k= (Ncrit t?)/(ir*D) (42) 

ACCURACY OF THE RESULTS 

Gene ral 

Table II shows the effect of grid size upon the error of the buck- 
ling coefficient for several cases. Table III shows the effect of as- 
pect ratio upon error for a single case. In many cases, the error is 
less than one per cent for rather coarse grid sizes. 

Factors Affecting Error 

General observations concerning the results are as follows: 

1. Error is greater for coarser grids, as expected. 

2. Error is greater for pure shear loading. 

3. Error is greater for clamped edges. 

4. Error is greater for narrower or longer plates. 
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TABLE II 



COMPARISON OF BUCKLING COEFFICIENTS WITH KNOWN VALUES FOR 
SEVERAL CASES, USING SEVERAL GRID SIZES 



Uniaxial Compression 









Aspect Ratio 


= 1 










Poisson's Ratio = 0 






ss 




k(known) = 4. 


00000 
















ss 




Qc 


Grid Size 


k 


7o Difference 






J J 


2x2 


4.01575 


+0.4 








3x3 


4.00324 


+0.08 




ss 




4x4 


4.00104 


+0.03 








6x6 


4.00021 


+0.005 








8x8 


4.00007 


+0.002 








Isotropic Compression 






I Ll i 




Aspect Ratio 


= 1 






4 4SS4 t 




PrW c e on 1 s Pot”} n = 0 










k(known) = 2. 


00000 
















SS 




ss 


Grid Size 


k 


7o Difference 




















2x2 


2.00791 


+0.4 








4x4 


2.00052 


40.03 








8x8 


2.00003 


+0.002 








Isotropic Compression 










Aspect Ratio 


= 1 






*1 In i 1 




Poisson's Ratio = 0 








■ 


5.30 < k (known) <C 5.33- 


-lower 








bound used for comparison. 


a 




CL 














Grid Size 


k 


7. Difference 
















| *CLf 1 




2x2 


5.36223 


+1.2 




II II 




3x3 


5.36599 


+1.2 








4x4 


5.32713 


40.5 








6x6 


5.30900 


40.2 








8x8 


5.30544 


40.1 
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TABLE II (continued) 




SS 



Pure Shear 
Aspect Ratio = 1 
Poisson's Ratio = 0 
k (known) = 9.34 

Grid Size k % Difference 



2x2 10.01581 +7.2 
3x3 9.57700 +2.5 
4x4 9.41822 +0.8 




Uniaxial Compression 
Aspect Ratio = 1 
Poisson's Ratio = 0.25 
k(known) = 1.70 

Grid Size k % Difference 



2x2 1.71399 +0.8 

4x4 1.69942 * 

8x4 1.69895 * 



*Zero for precision of known value. 




Pure Shear 
Aspect Ratio = 1 
Poisson's Ratio = 0 
k(known) = 14.71 

Grid Size k % Difference 



2x2 

3x3 

4x4 



23.26396 

16.04630 

15.04271 



+58.2 

+9.1 

+2.3 
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TABLE II (Continued) 




Uniaxial Compression 
Aspect Ratio = 1 
Poisson's Ratio ■ 0.25 
k (known) = 1.43418 

Grid Size k 



2x2 1.44324 
4x4 1.43483 
8x4 1.43435 



% Difference 

+ 0.6 

+0.05 

+ 0.01 



— CL— — 




Pure Shear 
Aspect Ratio = 1 
Poisson's Ratio = 0 
k(known) = 12.28 

Grid Size k 7» Difference 

2x2 17.13879 +39.6 

4x4 12.82566 +4.4 




Uniaxial Compression 
Aspect Ratio = 1 
Poisson's Ratio = 0 
k (known) = 6.74 

Grid Size k 

2x2 6.82334 

4x4 6.78242 

8x8 6.74583 



7« Difference 

+ 1.2 
+0 . 6 
+0.09 





CL 




Uniaxial Compression 
Aspect Ratio = 1 
Poisson's Ratio = 0 
k (known) = 7.69 




js 




ss 


Grid Size k 


7o Difference 








2x2 8.79292 


+14.3 




CL 




4x4 7.74665 

8x8 7.69541 


+0.7 

+0.07 
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TABLE III 



COMPARISON OF BUCKLING COEFFICIENTS WITH KNOWN VALUES 
FOR ONE CASE, USING SEVERAL ASPECT RATIOS 








Uniaxial Compression 
Grid Size, 4x4 
Poisson’s Ratio = 0 
Aspect ratios for buckling mode 
transition: 2% 3^, 4^, etc. 



a/b 


k 


k (known) 


Z Error 


No. Half Waves 


0.2 


27.05296 


27.04000 


+0.048 


1 


0.4 


8.41330 


8.41000 


+0.039 


1 


0.6 


5.13940 


5.13778 


+0.032 


1 


0.8 


4.20364 


4.20250 


+0.027 


1 


1.0 


4.00104 


4.00000 


+0.026 


1 


1.2 


4.13556 


4.13434 


+0.030 


1 


1.4 


4.47150 


4.47022 


+0.029 


1 


2.0 


4.00770 


4.00000 


+0.193 


2 


3.0 


4.02924 


4.00000 


+0.731 


3 
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When the buckled shapes of the plates are considered, all of these 
observations may be further generalized to a single one. The accuracy 
of representing the buckled surface is less for greater curvature gradi- 
ents. To keep error low, finer grids must be used for more complex 
buckled shapes. 

Sources of Error 

Discretization Error Most numerical methods replace continuous 
equations with discrete approximations. The finite element method, 
however, uses continuous equations but a discrete approximation of the 
continuous structure being analyzed. The finite element method has in- 
herent discretization error. Since the finite element approximation 
of a real plate is stiffer than the real plate, the discretization 
error of the buckling coefficient is never negative. In addition, the 
discretization error decreases montonically as the plate is further 

subdivided, provided that previously existing nodes are also nodes in 
2 5 

the finer grid. * For errors shown in Tables II, III, and IV, most 
error is believed to be discretization error. 

Manipulation Error Matrix iteration produced monotonically con- 
verging values of the buckling coefficient. Further, as a consequence 
of the matrix iteration method used, convergence was always from an 
initially higher value. Hence, matrix iteration error of the buckling 
coefficient is also positive. A convergence criterion of .001 per cent 
difference between successive values of the buckling coefficient was 
used. Consequently, the error of the results cannot be much less than 
.001 per cent. Round-off errors are believed to be insignificant in 
comparison with other errors, since double precision programming was 
used . 
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Negative Error 



Because the two principal sources of error--discretization and 
convergence --both produce positive error of the buckling coefficient, 
the error cannot be negative. Although some errors initially appeared 
to be negative , in each case encountered the cause was found to be lack 
of precision of the "known 1 ' value of the buckling coefficient. 

Accuracy of Buckling Mode Vectors 

Although the buckled shape of the plate is of secondary interest, 
the buckling mode vector was available as a consequence of the method used. 
The vector provided another check on the accuracy of the method. For 
example, the nodal displacements were compared with values computed from 
the known buckled shape for the uniaxially compressed, simply supported 
plate of aspect ratio one for a 4 x 4 grid. The displacements were in- 
ternally consistent to a minimum of three significant figures. 

EFFECT OF SIXTEEN-DEGREE -OF -FREEDOM ELEMENT 

Table IV shows comparisons of results with those given by Kapur and 
4 - 

Hartz for two cases. An improvement in accuracy is apparent. More 
important, for the same order of magnitude of the error, fewer elements 
were required. 

To show what this improvement means, consider the case of the simply 
supported plate compressed in one direction. Kapur and Hartz report an 
error of -0.58 per cent when using a 12 x 12 grid. Assuming that no 
symmetry and no compaction of matrices is used, and not considering re- 
duction of the matrix order when boundary conditions are applied, their 
formulation requires a 507 x 507 assembled stiffness matrix to represent 
three degrees of freedom at each of 169 nodes. A stability coefficient 
matrix of the same order is also required, of course. Using one more 



30 



TABLE IV 



COMPARISON OF BUCKLING COEFFICIENTS WITH VALUES OBTAINED 
USING TWELVE -DEGREE-OF-FREEDOM ELEMENT 3 





SS 






Uniaxial Compression 








— 




Aspect Ratio = 1 












Poisson 1 s 


Ratio = 0 




ss 




SS 




k (known) : 


= 4.00000 
























12 Degrees 


of Freedom 


16 Degrees 


of Freedom 


Grid Size 




k 


% Error 


k 


7. Error 




2x2 




***** 


**** 


4.01575 


+0.39 




3x3 




3.645 


-8.88 


4.00324 


40.08 




4x4 




3.77 


-5.75 


4.00104 


40.03 




6x6 




3.887 


-2.80 


4.00021 


40.005 




8x8 




3.933 


-1.68 


4.00007 


40.002 




10 x 10 




3.96 


-1.00 


******* 


****** 




12 x 12 




3.977 


-0.58 


******* 


****** 




1 ki-j J 






Isotropic 


Compression 












Aspect Ratio = 1 





Poisson's Ratio = 0 
5.30< k (known )< 5.33 C 





-* — 12 Degrees 


of Freedom 


16 Degrees 


of Freedom 


4 f f 

Grid Size 


k 


1 Diff. 


k 


7. Diff 


2x2 


4.901 


-7.53 


5.36223 


+1.17 


3x3 


***** 


***** 


5.36599 


+1.24 


4x4 


4.975 


-6.13 


5.32713 


40.51 


6x6 


5.078 


-4.19 


5.30900 


40.17 


8x8 


5.160 


-2.64 


5.30544 


40.10 


10 x 10 


5.217 


-1.57 


******* 


***** 



Data from Kapur and Hartz, reference 4. 

^The lower bound was used for computing per cent differences. 

Since Kapur and Hartz had used an average of the upper and lower bounds, 
the per cent differences given by them were re-computed. 
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degree of freedom per node, the error is 0.39 per cent when using a 
2x2 grid. For the same assumptions, the stiffness matrix size is 36 x 
36 for four degrees of freedom at each of nine nodes. In this example, 
the superior accuracy of the newer plate element and the consequent 
savings in computer time and space are unequivocally apparent. 

PREVIOUSLY UNSOLVED PROBLEMS 

Table V shows the results of investigation of two previously un- 
solved plate buckling problems. 

An edge which has more than one edge condition along its length 
is difficult to analyze by other means. Although the computer program 
used for other cases had to be altered very slightly to handle such an 
edge, there was no particular difficulty in using the finite element 
method. 
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TABLE V 



PLATE BUCKLING PROBLEMS NOT PREVIOUSLY INVESTIGATED 







FR 


FR 


X 


X 


X 














1 


T \ 




^4 


5' 
















i 


1 


1 

| 














FR 


<C> 


1 

jT 


T® 


1 

Is 


lO 


•X 














1 


i 


1 














■>* 

FR' 


11 


i 

\iZ 

\ 


\ 

il5 


1 

*14 

f 


45 


X 






















F S' 




" ~T lT " 

i 


*18 


[IS 

!z4 


20 


X 












1 

1 25 


25 














21 


122 












F 


R 


FR 


X 


X 




X 






Uniaxial compression, 


aspect 


ratio = 2, 


Poisson's ratio = 0.3 






Nodes 


marked 


X are either all 


simply supported or all clamped. 












HALF FREE, 


HALF SIMPLY SUPPORTED 




















k = 0.179 
















Relative 


Transverse Displacements of Nodes 






Node 


w 




Node 


w 


Node 


w Node 




w 


Node 


w 


i 


-.3562 




2 


-.1207 


3 


0 4 




0 


5 


0 


6 


-.3662 




7 


-.1354 


8 


.0139 9 


+.0019 


10 


0 


11 


-.3699 




12 


- . 1408 


13 


.0182 14 


+.0027 


15 


0 


16 


-.3662 




17 


-.1354 


18 


.0139 19 


+.0019 


20 


0 


21 


-.3562 




22 


-.1207 


23 


0 24 




0 


25 


0 



Node 


w 


HALF FREE, HALF CLAMPED 
k = 0.208 

Relative Transverse Displacements of 
Node w Node w Node 


Nodes 

w 


Node 


w 


i 


-.3174 


2 


-.0956 


3 


0 


4 


0 


5 


0 


6 


-.3257 


7 


-.1073 


8 


-.0034 


9 


+.0013 


10 


0 


11 


-.3287 


12 


-.1111 


13 


-.0053 


14 


+.0026 


15 


0 


16 


-.3257 


17 


-.1073 


18 


-.0034 


19 


+. 0013 


20 


0 


21 


-.3174 


22 


-.0956 


23 


0 


24 


0 


25 


0 
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4. THE COMPUTER PROGRAM 



GENERAL 

A sample computer run is given in Appendix I. The program listing 
is largely self-explanatory. The program is in FORTRAN IV with double 
precision. The IBM 360 computer at the Naval Postgraduate School was 
used. This computer has a core memory capacity of 512,000 bytes. The 
HASP compiler used required about 55,000 bytes. 

CAUTIONS 

Program Efficiency 

The program was written as an aid in producing this thesis only, 
by a relatively inexperienced programmer. Paradoxically, one purpose 
of the thesis is to show how computer time and space can be saved, yet 
both were often sacrificed for ease of programming. 

Convergence in Matrix Iteration 

Initially, a trial vector with a high content of some mode shapes 
was used to start the matrix iteration process. This sometimes caused 
convergence to an incorrect mode. So that the trial vector would not 
be prejudiced toward particular modes, a new trial vector consisting of 
elements taken from a table of random numbers was used. In all cases 
run using the new trial vector, where the proper buckling mode was 
known, convergence was to the correct mode. 

Convergence was very slow for cases with the two largest eigen- 
values nearly equal. Although methods of accelerating convergence were 
tried, no generally acceptable method was found. The program is not 
recommended for such cases. 
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5. SUMMARY 



CONCLUSIONS 

For plate buckling problems using the finite element method, 
the sixteen-degree-of -freedom element is far more accurate than the 
twelve-degree-of -freedom element. Not only is the accuracy superior, 
but it is sufficiently greater so that sizes of assembled stiffness 
and stability coefficient matrices may be greatly reduced. Far fewer 
plate elements are needed for equivalent accuracy. Computer time and 
space requirements are significantly reduced. 



PROBLEMS NEEDING FURTHER STUDY 

The following problems were beyond the scope of this investigation: 

1. An attempt should be made to find a general rule for sub- 

dividing the plate for minimum error. Elements of un- 
equal sizes, or more divisions in one direction might be 
used, for example. 

2. A systematic investigation of many of the remaining unsolved 

cases of plate buckling should be made. 

3. Further study of methods of determining the eigenvalue 

should be made in order to find a more rapid method. 

4. A more general and more efficient computer program 

should be devised for industrial use. 
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APPENDIX I 
SAMPLE COMPUTER RUN 
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non noon noon non oooooooooono oonooo 



MS ^V!lfi?I A gLl»E» n mB 0 Bi?H c Si L i8 B x 0 S6 R i£I«g8V L «?iPWiy» 

MATRIX. 



MAXIMUM PLATE SIZE A ELEMENTS X 4 ELEMENTS OR EQUIVALENT. 



1 IMPLICIT REAL * 8 (A-H,0-Z) 

2 DIMENSION POL Y 1 ( 16 > , POL Y2 ( 16 ) f POL Y3 (1 6 ) , POL Y4 ( 1 <S ) . POL Y5 ( 1 6 ) f 
1P0LY6( 16) f A( 16,16) ,AINV( 16, 16) ,SELMT( 16,16),P(3,3), GCOFF( 3,16), 
2GX EXP ( 3, 16) • GYEXP ( 3 , 16 ) , G I NT ( 1 6. 16 ) . S ( 1 00 . 1 00 ) , $ I NV ( 100 , 1 00 ) , 
3BFLMT( 16,16) • HC0EF ( 2 , 1 6 ) , HXE XP ( 2 . 1 6 I . HYEX P ( 2 , 1 6 ) , H I NT ( 1 6 . 1 6 ) • 

48 ( 100, 100) ,C( 2.2) ,W( 100) ,WNE XT ( 100) , STIFF ( 100, 100) ,WSAVE( 100) 

201 EQUIVALENCE (A(1),GINT(1),HINT(1),B(1),S(1)) , ( GCOFF ( 1 ) • HCOEF ( 1 ) 
1 SELMT ( 1 ) • BEL MT ( 1 ) , WSA VE ( 1) ) , ( GXE XP ( 1) ,HXEXP( 1 ) ) « ( GYEXP I 1) « 
2HYFXPI 1 ) > , ( W( 1 ) ,P0LY1( 1 ) > , ( WNE XT < 1 ) , POLY2 ( l ) ) , ( P (1) , ST I FF ( 1 ) ) 



INPUT FROM FIRST DATA CARD. ASPFCT RATIO (AR), POISSON'S RATIO 
(POISS), NO. ELEMENTS IN X DIRECTION ( MEL MT ) , NO. ELEMENTS IN Y 
DIRECTION (NELMT), BOUNDARY CONDITIONS FOR Y = 0 , Y = 8, X = 0, X = A 
( IBND,JBNO,KBND,LBND) , RELATIVE SIZES OF LOADS NX, NY, NX Y (Cl, 
C2 , C 3 ) • FORMAT 2D10. 0,613, 3D10.0. 



FOR BOUNDARY CONDITIONS, 1=FRE E, 2 = SIMPLY SUPPORTED, 3 = CL AMP E D , 
4= SYMMF TR Y , 5= ANTI -S YMMF TRY . 4 AND 5 MAY BE USEO ONLY FOR JBND 
AND LBND. NX AND NY MUST NOT BE NEGATIVF ( COMP R E SS I ON ONL Y ) . NX 
AND NY MUST BE ZERO IF NXY IS OTHER THAN ZERO (PURE SHEAR ONLY). 



196 READ (5.197) AR , PO I SS , MEL MT , NEL MT , I BNO, J BND , KBNO , L BND, C 1 , C2 , C 3 

197 FORMAT (2D10.0 ,613,3010.0) 

WRITE INPUT DATA 



198 WRITE (6,199) AR , P 0 I SS , MFL MT , NFL MT 

199 FORMAT { 1H 1 , 1 3H A SP ECT R A T I 0= , F 8 . 5 , 3X , 1 4H P 0 1 S SON R AT I 0= , F4 . 2 , 3X , 
11PHN0. FLE^FNTS IN X = , I 1 , 3 X » 1 8HNQ. ELEMENT S IN Y=,I1) 

202 WRITE (6.203) 

203 FORMAT ( lH 0 , 69HB0UNDAR Y CONDI T I ONS , 1 = FR E E , 2=PINNED, 3=CLAMPF0, 4 
1=SYM. , 5 = A NT I - SYM. ) 

WRITE (6,204) I BND , JBND , K BND , L BND 

204 FORMAT ( 1H 0 , 18HY=0 , EDGE COND I T ION, I 2, 3X , 1 8HY = B, EDGE CONDI T ION, I 2, 

1 3X, 1 8HX = 0 * EDGE CONDITION, I2,3X,1RHX=A,EDGE COND I T I ON , I 2 ) 

WRITE (6.3) C1,C2,C3 

3 FORMAT ( 1H0,24HRELATI VE LOAD SIZES NX= , FS . 2 , 2 X , 3HNY=, F5 . 2 , 2X , 
14HNXY=,F5. 2) 

THIS PORTION OF THE PROGRAM PRODUCES THE RECTANGULAR PLATE ELEMENT 
STIFFNESS MATRIX. 

205 EM=MELMT 

206 EN=NEL MT 
EX= 1 ./FM 
EY=l./( AR*EN) 

IF ( (LBND.EQ.4) .OR. (LBND.EQ.5) ) F X= 1 . / ( 2 . * F M ) 

IF ( (JBND.EQ.4).0R.( JBND.E0.5) ) EY=1 . / ( 2. *EN* AR ) 

ESTABLISH ARRAYS FOR F XPONENT S AND MULTIPLIERS FOR TERMS OF 
POLYNOMIAL W AND ITS DERIVATIVES 



4 K=1 

5 DO 13 1=1,4 

6 DO 13 J=l,4 

7 POL Y 1 ( K ) = I — 1 

8 POL Y2( K )= J- 1 

9 P0LY3(K ) = I - 2 

10 POL Y4( K ) = J -2 

11 POL Y 5 ( K ) = I - 3 

12 POLY 6 ( K ) = J -3 

IF ( POLYK K) .LT.O. ) 
IF ( POL Y2(K). LT.O.) 
IF (P0LY3(K) .LT.O. ) 
IF ( P0LY4( K) .LT.O. ) 
IF ( POLY 5 ( K ) .LT.O. ) 
IF ( POL Y6 ( K ) .LT.O. ) 

13 K=K*1 



P0LY1 ( K ) =0 . 
P0LY2(K)=0. 
P0LY3( K ) =0 • 
P0LY4 ( K ) =0 • 
POL Y5 ( K) =0 • 
POL Y6 ( K) =0. 



ESTABLISH TRANSFORMATION MATRIX A. 



15 DO 28 J = 1 , 1 6 

16 X=0 . 

17 IF ( ( I.GT.4.AND. I.LT.9).0R. ( I.GT.12) ) X=FX 

18 Y= 0 . 

19 IF (I . G T . 8 ) Y=E Y 

20 IF ( ( X.FO.O. ) .AND. ( P0LY1 ( J ) .EQ.O. ) ) GO TO 23 

21 TERM1=X**P0LY1 ( J ) 

22 GO TO 24 

23 TERM 1= 1 . 

24 IF ( ( Y.EO.O. ) .AND. ( P0LY2( J) .EO.O. ) ) GO TO 27 

25 TERM2=Y**POLY2(J ) 

26 GO TO 28 

27 TERM2= 1 . 

28 A ( I , J )=TERM1*tERM2 

29 DO 43 1=2,14,4 

30 DO 43 J=l, 16 

31 X=0 • 

32 IF ( ( I .GT.4. AND. I .LT.9) .OR. ( I . GT . 12) ) X=EX 

33 Y= 0 . 

34 IF ( I .GT.8 ) Y=EY 

35 IF (( X. EO.O. ) .AND. ( P0LY1 ( J) . EQ.O.) ) GO TO 38 

36 TERM1=X**P0LY1(J) 

37 GO TO 39 

38 T FRM 1 = 1 . 

39 IF ( (Y. EQ.O. ). AND. ( POL Y4( J) .EQ.O. ) ) GO TO 42 

40 TERM2=Y**P0LY4(J ) 

41 GO TO 43 

42 TERM2=1. 

4} A(I,J)=P0LY2(J)*TERM1*TFRM2*EY 
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44 00 58 1=3,15,4 

45 00 58 J = 1 , 16 

46 X= 0 • 

47 IF ( ( I.GT.4.AND. l.LT.9).0R.( I.GT.12) ) X=FX 

48 Y= 0 . 

49 IF ( I .GT.R ) Y=E Y 

50 IF ( ( X.E0.0. ) . AND. ( POL Y3( J) .EQ.0. ) ) GO TO 53 

51 TERM1=X**PDLY3( J) 

52 GO TO S4 

53 TERM 1= 1 . 

54 IF ( (Y.E0.0.).AN0.(P0tY2( JJ.EG.O.)) GO TO 57 

55 TFRM2=Y**P01Y2< J) 

56 GO TO 5 8 

57 TERM2= 1 • 

58 A(I,J)=P0LY1(J)*TERM1*TERM2*EX 

59 00 73 1=4, 16,4 

60 00 73 J = l, 16 

61 X= 0 • 

62 IF ( ( I .GT.4.AN0. I .IT.9) .0R.( I.GT.12 ) ) X=FX 

63 Y= 0 • 

64 IF ( I.GT.8) Y=FY 

65 IF ( ( X.E0.0. ) .AND. ( P0LY3( J J.E0.0. ) ) GO TO 68 

66 TERM1=X**P0LY3(J) 

67 GO TO 69 

68 TFRM 1= 1 . 

69 IF ( ( Y.EO.O. ) .AND. (P0LY4( J ) . EQ.O. ) ) GO TO 72 

70 TERM2=Y**POLY4( J) 

71 GO TO 73 

72 TERM2=1. 

73 A ( I V J)«P0LY1(J)*P0I.Y2( J)*TERM1*TERM2*EX*EY 

FIND INVERSE OF MATRIX A (AINV). GAUSS-JORDAN METHOD. 

86 00 90 1=1,16 

87 DO 90 J = 1 , 16 

88 A l N V ( I, J)=0. 

89 IF ( I .EQ.J) A I NV ( I , J) = l. 

90 CONTINUE 

91 EP S= 1 .0- 10 

92 DO 113 1=1,16 
9 3 K= I 

94 IF (1-16) 95,102,95 

95 IF ( A( I , I ) -EPS ) 96,97, 102 

96 IF ( -A ( 1,1 )-FPS) 97,97,102 

97 K= K+ 1 

98 DO 100 J= 1,16 

99 AINV( I, J)=AINV( I , J ) +A I NV ( K , J ) 

100 A ( l ,J)=A( I , J ) ♦ A ( K , J ) 

101 GO TO 95 

102 D I V= A ( 1,1) 

103 DO 105 J= 1 , 1 6 

104 AINV( I , J)=AINV( I , J)/DI V 

105 A ( I , J ) = A ( I , J )/DI V 

106 DO 113 1=1.16 

107 DELT = Ml. I ) 

108 IF (DABS(DELT)-EPS) 113,113,109 

109 IF (L-I ) 110,113,110 

110 DO 112 J= 1 , 16 

111 AINV(L,J)=AINV(L.J)-AINV(I , J)*DELT 

112 A(L,J)=A(L,J)-A( I,J)*DELT 

113 CONTINUE 

DO 843 1=1,16 
DO 843 J= 1 , 1 6 

IF (DABS( UNV( I, J) ).LT. 1.0-10) AINV(I,J)=0. 

843 CONTINUE 

ESTABLISH MATRIX P, THE MATRIX OF FUNCTIONS OF POISSON'S RATIO 



126 P( 1, 1) = 1. 

127 P( 1 , 2 ) =P0 1 SS 

128 P< 1 , 3 ) = 0 . 

129 P ( 2 , 1)=P0ISS 

130 P ( 2 , 2 ) = 1 . 

131 P< 2 , 3 ) = 0 . 

132 P( 3, 1)=0. 

133 P( 3, 2)=0. 

P(3,3)=2.*( l.-POISS) 

ESTABLISH MATRIX G 



134 DO 143 1=1 

135 GCOE F ( 1 , I ) 

136 GCOE F ( 2 , I ) 

137 GCOE F ( 3 , I ) 

138 G X E X P ( 1, I ) 

139 GXEXP( 2, I ) 

140 GXEXP ( 3 , I ) 

141 GYEXP ( 1,1) 

142 GY F XP ( 2 » I ) 

143 GYEXP (3,1 ) 



,16 

=PPLY1( I )*P0LY3( I) 
=PCLY2( I ) *P0LY4( I ) 
= P0LY1 ( I ) ♦POL Y 2 ( I ) 
=P0LY5( I ) 

= PCL Y 1 ( I ) 

= PCLY3( I ) 

=PPLY2( I ) 

= P0LY6( I ) 

= P0LY4( I ) 



THE FOLLOWING PART OF THE PROGRAM COMPUTES I NT FGR \[ OVFP X FRhm 
0 TO EX, INTFGRAL OVER Y FROM 0 TO FY, OF G TRANSPOSE ♦ P ♦ G. THF 
RESULT IS NAMED GINT. 

144 DO 162 1=1,16 

145 DO 162 J= 1,16 

146 SUM1=0. 

147 DO 151 K= 1,3 

148 SUM2=0. 

149 DO 150 L= 1,3 

150 SUM2=SUM2*GCOEF( K, I ) *P( K,l ) ♦GCOE F ( L , J ) * ( 1 . / ( ( GXFXP ( K . I ) + GXEXP < L , J ) 
1 + 1 . ) ♦ ( GYF XP ( K , I)+GYEXP(LtJ) + l. ) ) )♦( E Y ♦♦ ( GY P X P ( K f I ) +GY FXP ( L , J ) + 1 . ) ) 
2*(EX**(GXEXP(K«I ) +GXE XP(L,J) + 1.) ) 

151 SUM1=SUM1+SUM2 
162 G I NT ( I, J)=SUM1 
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FIND ELEMENT STIFFNESS MATRIX (SELMT). SfrL MT= A I N V TR ANSP0SF*G I NT* 
A I N V 

175 00 183 1=1,16 

176 00 183 J= 1,16 

\n aH M ti2-K=i,i6 

179 SUM2=0. 

180 00 181 L = 1*16 

181 $UM2=SUM2*GINT(K,L ) ♦A I NV( l f J) 

182 SUM 1 = SUM 1 ♦ A I NV ( K , I ) *SUM2 

183 SE LMT ( I , J ) =SUM 1 

GENERATE ASSEMBLED STIFFNESS MATRIX S. THE FIRST ELEMENTS OF 4 X 4 
ARRAYS ARE ASSIGNED THEIR PROPER LOCATIONS. THE REMAINDER OF EACH 
ARRAY IS GENERATED FROM THE FIRST ELEMENT. 

200 INDEX=( MELMT + 1 ) *( NELMT + 1 ) *4 
I 1 = 1 

209 DP 211 1=1, INDEX 

210 DO 211 J= 1 , 1 NDEX 

211 S ( I , J ) = 0 . 

212 M=M§LMT+1 
N=NELMT + 1 

214 11=0 

215 DO 372 I = 1 » M 

216 DP 372 J*1,N 

217 J J=0 

218 DO 371 K= 1 » M 

219 DO 371 L = 1 , N 

2 20 IF ( (K.GT. ( IU) ) .OR. (K.LT* ( 1-1 ) ) .0R.( L.GT. ( J + l ) ) .OR. ( L.LT. ( J-l) ) ) 

1 GO TO 371 

221 IF (K-I) 222,227,228 

222 IF (L-J) 233,238.245 

2 33 IF ( ( I .EQ. 1) .OR. ( J.EO. 1) ) GO TO 371 
GO TO 224 

238 IF ( I . EO . 1 ) GO TO 371 
IF (J.EO.l) GO TO 230 
IF (J.EO.N) GO TO 235 
GO TO 240 

245 IF <<I .EG. 1) .OR. < J .EO.N) ) GO TO 371 
GO TO 248 

227 IF (L-J) 246,251.252 

246 IF (J.EO.l) GO TO 371 
IF ( I .EO. 1 ) GO TO 254 
IF ( I.EO.M) GO TO 259 
GO TO 264 

251 IF ( I.EO. 1) GO TO 257 
IF (I . EO. M ) GO TO 262 
IF (J.EO.N) GO TO 298 
IF (J.EO.l) GO TO 305 
GO TO 319 

262 IF (J.EO.N) GO TO 288 

IF (J.EO.l) GO TO 293 

GO TO 312 

257 IF (J.EO.l) GO TO 271 

IF (J.EO.N) GO TO 276 

GO TO 281 

252 IF (J.EO.N) GO TO 371 
IF ( I .EO. 1 ) GO TO 394 
IF (I.EO.M) GO TO 397 
GO TO 400 

228 IF (L-J) 269,274,279 

269 IF ( ( I .EO.M) .0R.( J.EO. 1) ) GO TO 371 
GO TO 405 

274 IF ( I .EO.M) GO TO 371 

IF (J.EO.l) GO TO 408 

IF (J.EO.N) GO TO 411 

GO TO 414 

279 IF (( I .EO.M) .OR. (J.EO.N) ) GO TO 371 
GP TO 419 

224 111=13 

225 JJJ=1 

226 GO TO 364 

230 I I I = 5 

231 JJJ=1 

232 GO TO 364 

235 111=13 

(lrf J TO 364 

240 1 1 1 = 5 

241 J J J = 1 

242 KK K= 1 3 

243 LLL=9 

244 GO TO 354 

248 I I 1 = 5 

249 J J J= 9 

250 GO TO 364 

254 I I 1 = 9 

255 JJJ=1 

256 GO TO 364 

259 111=13 

260 J J J= 5 

261 GO TO 364 

264 I I 1 = 9 

265 JJJ=1 

266 KK K= 1 3 

267 LLL=5 

268 GO TO 354 

271 I I 1 = 1 

272 JJJ=1 

273 GO TO 364 

276 I I 1 = 9 

277 J J J = 9 

278 GO TO 364 

281 I I I = 1 

282 JJJ=1 

283 KK K= 9 
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284 


LLL = 9 




285 


GO TO 


354 


288 


m = i3 




289 


J J J* 1 3 




290 


GO TO 


364 


293 


I I 1*5 




294 


J J J = 5 




295 


GO TO 


364 


298 


I I 1 = 9 




299 


J J J = 9 




300 


K KK= 1 3 




301 


LLL* 1 3 




302 


GO TO 


354 


305 


I I 1*1 




306 


JJJ=1 




307 


KK K-= 5 




308 


LLL=5 




309 


GO TO 


354 


312 


I I 1*5 




313 


JJ J=5 




314 


KK K= 1 3 




315 


LLL* 1 3 




316 


GO TO 


354 


319 


I 11 = 1 




320 


JJJ*1 




321 


K K K= 5 




322 


LL L=5 




323 


MMM* 9 




324 


NNN* 9 




325 


I I I 1=13 


326 


JJ JJ=13 


327 


GO TO 


340 


394 


I I I = 1 




395 


J J J *9 




396 


GO TO 


364 


397 


I 11 = 5 




398 


JJ J=13 




399 


GO TO 


364 


400 


l I 1=1 




401 


JJ J=9 




402 


KK K= 5 




403 


LL L= 1 3 




404 


GO TO 


354 


« l JJJ:? 




407 


GO TO 


364 


408 


I I 1 = 1 




409 


JJ J = 5 




410 


GO TO 


364 


411 


I I 1*9 




412 


JJ J=13 




413 


GO TO 


364 


414 


I I 1=1 




415 


JJ J=5 




416 


KK K = 9 




417 


LLL* l 3 




418 


GO TO 


354 


419 


I I 1 = 1 




420 


JJ J=13 




421 


GO TO 


364 


340 


DO 352 


KK 


341 


MM=I I*KK 


342 


DO 348 


LL 



343 NN = J J +LL 

344 SI MM,NN)=SELMTl I I I , JJJ J+SFLMTI KKK, LLL )+SELMT<MMM f NMN) 
1 +SELMT ( I I I I * J J J J ) 

345 JJJ=JJJ+1 

346 LLL=LLL+ 1 

347 NNN=NNN+ 1 

348 JJJJ=JJJJ+1 
J J J* J J J-4 
LLL=LLL-4 
NNN=NNN-4 
JJ JJ=JJJJ-4 

349 I I 1 = 1 I I + l 

350 KK K=KKK 4- 1 

351 MMM=MMM+1 

352 II 11=1111+1 

353 GO TO 371 

354 DO 362 K K= 1 , 4 

355 MM = I I +K K 

356 DO 360 LL=1,4 

357 NN=JJ+LL 

3 58 S(MM t NN) = SELMT( I II , J J J ) + SE LMT ( KK K t L ID 

359 JJJ=JJJ+1 

360 LLL=LLL + 1 
J J J = J J J-4 
LLL=LLL-4 

361 111=111+1 

362 KK K = KKK 1 

363 GO TO 371 

364 DO 370 KK=1,4 

365 MM= I I + KK 

366 DO 369 LL*1,4 

367 NN= J J+Ll 

368 SI MM,NN)=SFLMTI I I I , J J J ) 

369 JJJ=JJJ+1 
J J J=J J J-4 

370 I I 1=1 I IM 

371 JJ=JJ+4 

372 11=11+4 

DO 844 1*1, INDEX 
DO 844 J* 1 , I NDEX 

IF (DABSISI l,J)) ,LT. 1.0-10) SII,J)=0. 

844 CONTINUE 

USE BOUNDARY CONDITIONS TO MOD I F Y S. 
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527 IF UIBNO.GT. 2). OR.UBNO.LT. 2)1 GO TO 529 
00 528 1*1, M 

n=4*(l-l)*N + l 
JJ-II+2 

DO 528 KK= 1.1 NDE X 
S< II,KK)*0. 

S(KK,II)*0. 

S( JJ.KKI-O. 

U K »HV!I8 ^KK) .AND. < Il.EO.l) ) S( I1,KK) = 1. 
IF ( ( JJ.EQ.KK) .AND* ( II .60,1) ) $(JJ,KK)=l. 

528 CONTINUE 

529 IF ( I BND • L T . 3) GO TO 531 
00 530 I = l , M 

1 i = 4 * ( i-n*N*i 
JJ=l 1 + 3 

00 530 KK=1, INDEX 
00 530 LL = I I » J J 
S ( LL ,KK ) = 0 • 

S(KK.LL)=0. 

IF < (LL. EQ.KK) .AND. (II .EQ. 1) ) S(LL,KK)=l. 

530 CONTINUE 

531 IF ( ( JBND.EQ.2).0R.(JBND.EQ.5)) GO TO 543 
GO TO 533 

543 00 532 1 = 1, M 
I l=4*I*N-3 
JJ=IU2 

00 532 KK= 1 , INOEX 
S( I I , KK ) = 0 • 

S ( KK , I I ) =0 • 

S( JJ,KK)=0. 

S(KK.JJ)=0. 

IF ( ( II .EQ.KK) .AND.( Il.EO.l) ) S(II,KK)=i. 
IF ( (JJ.EQ.KK) .ANO.C Il.EO.l)) S(JJ,KK)=l. 

532 CONTINUE 

533 IF < ( JBND.GT.3) .OR. ( JBN0.LT.3) ) GO TO 535 
00 534 1 = 1 , M 

I I=4*I*N-3 
JJ=lI+3 

00 534 KK=1, INDEX 
00 534 LL= I l , J J 
S ( LL , KK )=0 • 

S(KK.LL)=0. 

IF ( (LL.EQ.KK) . AND. ( 1 1 .EQ. 1) ) S(LL,KK) = 1. 

534 CONTINUE 

535 IF ( (JBNO.GT. 4). OR.UBNO.LT. 4) ) GO TO 537 
00 536 1=1, M 

1 1=4*1 *N- 2 
J J= 1 1 + 2 

00 536 KK= 1 f I NDE X 
S( l I , KK ) =0 • 

S ( KK , II)=0, 

S( JJ , KK 1=0. 

S(KK. JJ )=0. 

IF ( ( II .EQ.KK) .ANO.( 11 .EO.i) ) S(II,KK)=i. 

I F C C JJ.EQ.KK). A NO. ( I 1 . EO . 1) ) S( JJ,KK)=l. 

536 CONTINUE 

537 IF ( (KBND.GT.2).0R.(KBND.LT.2>) GO TO 539 

00 538 1 = 1, N 

1 1 = 4* 1-3 
JJ=I 1+1 

00 538 KK=1, INOEX 
S( II ,KK ) =0 • 

S( KK, I I ) =0 • 

S( JJ ,KK )=0 • 

S ( KK • J J ) =0 • 

IF <( II. EQ.KK) .AND. ( 11. EO. 1U S(II,KK)=l. 
IF ( (JJ.EQ.KK) .AND. ( Il.EO.l)) S(JJ,KK)=l. 

538 CONTINUE 

539 IF (KBN0.LT.3) GO TO 541 
00 540 1 = 1, N 

I I=4*I-3 
JJ=IU3 

00 540 KK=l f INDEX 
00 540 LL=I I , JJ 
S( LL,KK)=0. 

S( KK *LL ) = 0 • 

IF (ILL. EQ.KK). AND. (II. EO. l>) S(LL,KK) = l. 
CONTINUE 

IF ( (LBN0.E0.2).0R.(LBND.E0.5) ) GO TO 542 
GO TO 545 

W-UtHilUfU.l-3 

jj=ii+i 

00 544 KK= 1 , INOEX 
S( II,KK)=0. 

S ( KK , II 1 = 0 . 

S( JJ,KK)=0. 

S ( KK , J J ) =0 • 

IF ( ( I I .EQ.KK) .AND. ( II. EO. 1) ) S(II,KK) = l. 
IF (( JJ.EQ.KK) .AND. < II. EO. 1) ) S( JJ,KK)=1. 

544 CONTINUE 

545 IF ( (LBNO.GT. 3). OR.ILBNO.LT. 3) ) GO TO 547 
DO 546 1=1, N 

I I=4*( M-l ) *N+4*I-3 
J J= I I + 3 

DO 546 KK=1. INOEX 
00 546 LL= l I ,JJ 
S( LL , KK ) = 0 . 

S( KK.LL)=0. 

IF ( (LL. EQ.KK) .AND. ( Il.EO.l) ) S(LL,KK)=1. 

546 CONTINUE 

547 IF ( (LBN0.GT.4).0R. (LBND.LT.4) ) GO TO 548 
DO 548 1=1, N 

II=4*(M-l)*N+4*I-l 
JJ=I m 

DO 548 K K = 1 , INDEX 
S( I I , KK ) =0 • 



540 

541 

542 
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S I KK * I I )=0. 

SI JJ,KK>=0. 

SIKK. JJ)=0. 

IF ( III .FQ.KK) .AND. I II .EQ. 1) ) SIIl,KK)=l. 
IF! I JJ.EQ.KK) .AND. I I l.EQ. 1 ) ) S( JJ»KK)-1. 
548 CONTINUE 

IF ( I1.E0.2) GO TO 496 

SAVE S BY STORING IN ARRAY NAMED STIFF. 

DO 614 1 = 1, INDEX 
DO 614 J=1 , INDEX 
614 ST IFF ( I,J)=S!t,J) 

FIND S INVERSE (SINV). 



422 DO 426 1=1, INDEX 

423 DO 426 J=l, INDEX 

424 S I NV ( I * J ) =0 • 

425 IF ( I.EG.J) S I NV I I * J ) = 1 

426 CONTINUE 

427 DO 44 8 I = 1 , T NDEX 
42R K= I 

429 IF II-INDEX) 430,437,430 

430 IF (S(I.I)-EPS) 431,432,437 

431 IF < -S C I, I >-ERS) 432,432,437 

432 K 4- 1 

433 DO 435 J=l. INDEX 

434 SINV( I,J)=SINV(I,J)+SINV(K,J) 

435 S( I, J) = S( I ,J)*S(K,J) 

436 GO TO 430 

437 D I V= S( 1,1) 

438 DO 440 J=1,TNDEX 

439 SINVI I, J)=SINV(I ,J)/DIV 

440 S< I, J) = S( I ,J)/0IV 

441 DO 448 L=1*IN0EX 

442 DEL T = S ( L , I ) 

443 IF (DABSIDEI T)-EPS) 448,448,444 

444 IF (L-I) 445,443,445 

445 00 447 J=1.INDEX 

446 SINV(L,J)=SINV(L,J)-SINV(I ,J)*DELT 

447 S(L, J) = S(L, J ) - SC I ,J)*OELT 

448 CONTINUE 

DO 845 1=1, INDEX 
DO 845 J=l, INDEX 

IF ( DABS I SINVI I f J)).LT. 1.0-10) SINVI I , J) =0 . 

845 CONTINUE 

THIS PORTION OF THE PROGRAM PRODUCES THE RECTANGULAR o L ATF ELEMENT 
STABILITY MATRIX BY A PROCEDURE SIMILAR TO THAT FOR THF FLEMFNT 
STIFFNESS MATRIX 

ESTABLISH LOAD CONDITION MATRIX C. 



449 C< 1, 1) = C2 

450 C< 1 , 2 ) =C3 

451 C( 2, 1 ) = C 3 

452 C( 2 , 2 ) =C 1 



ESTABLISH MATRIX H 



453 DO 459 1=1,16 

454 HCOEFI 1,1 )=PCLY2< I ) 

455 HCOEFI 2,1 )=PCLY1 I I ) 

456 HXEXPI 1 , I )=POLYt I I ) 

467 HXEXPI 2,1 )=PCLY3( I ) 

458 HYEXPI 1,1 )=PCLY4( I ) 

459 HYEXPI 2,1 )=PCLY2< I ) 

THE FOLLOWING PART OF THF PROGRAM COMPUTES INTEGRAL OVER X FROM 0 
TO FX, INTEGRAL OVER Y FROM 0 TO EY, OF H TRANSPOSE * C * H. 

THE RESULT IS NAMED HINT. 

460 DO 468 1=1,16 

461 DO 468 J=1 ,16 

462 SUM1=0. 

463 DO 467 K=1 ,2 

464 SUM2=0. 

465 DO 466 L=1 ,2 

4 66 SUM2 = SUM2 + HCOEF( K, I ) *C I K , L ) *HC OFF ( L , J ) * I l./((HXFXP(K,I ) +HXFX P I L , J ) 
1*1 .) *1 HYEXPI K, I ) +HYFXP (L ,J)*1 .)))*< FY**I HYEXPI K, l ) *H Y EX P I I , J ) ♦ 1 . ) ) 
2*1 EX** ( HXEXPI K, I )*HXEXP(L,J) + 1. ) ) 

467 SUM1=SUM1*SUM2 

468 HINT! I , J ) = SUM1 

FIND ELEMENT STA8ILITY MATRIX IBELMT). BFLMT=AINV TRANSPOSE * 

HINT * AINV. 

476 DO 484 1 = 1 , !6 

477 DO 484 J=1 ,16 

478 SUM1 =0 • 

479 DO 483 K=1 , 16 

480 SUM2=0. 

481 DO 482 L= 1 ,16 

482 SUM2=SUM2*HINT(K,L )*AINV(L ,J) 

483 SUM1 = SUMUAINVIK,I )*SUM2 

484 BELMTI I , J )=SUMt 

GENERATE ASSEM8LED STA8ILITY MATRIX AND APPLY BOUNDARY CONDITIONS 
BY USING SAME PORTION OF PROGRAM USED FOR STIFFNESS MATRIX. 

494 11=2 

495 GO TO 209 

496 CONTINUE 



THIS PORTION FINDS THE RUCKLING COEFFICIENT AND OISPLACFMFNT 
VECTOR BY MATRIX ITERATION. 
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1 = 1 
MM= 1 

ESTABLISH TRIAL VECTOR OF ELEMENTS RETWEEN -1 AND *1 
TABLE OF RANODM OIGITS 

11 = 1 
JJ=10 

497 DO 739 1=1,10 

READ (5,740) (W( J) ,J=I l,JJ> 

11=11*10 

WRaE^(6l741 ) 2 INDEX 

741 FORMAT ( lH0,19HTRIAL VECTOR, F I R ST , I 4 , 1 X , 1 7HELEMENTS 
11 = 1 

JJ = 10 

DO 498 1=1.10 

WRITE (6,742) ( W( J ) , J=I I , J J) 

11=11*10 
JJ=J J* 10 

498 CONTINUE 

742 FORMAT (10F6.2) 

MULTIPLY W TRANSPOSE * S * W = WTSW. 

700 WTSW=0 • 

DO 701 1=1, INDEX 
00 701 J=l, INDEX 

701 WTSW=WTSW*W( I ) ♦ST I FF ( I , J ) *W( J ) 

NORMALIZE W. 

DO 702 1=1, INDEX 

702 W(I)=W( I)ZDSORT(WTSW) 

IF (MM.E0.1) GO TO 713 

AFTER FINAL I TER AT I ON, WR I TE TRANSVERSE DISPLACEMENTS 

722 IF (C3.EQ.0.) GO TO 721 
DO 716 1=1, INDEX 
IF (MM.EQ.2) W( I )=WSAVE( I )*W( I ) 

IF (MM. EG. 3) W(I )=2.*WSAVE( I )-W( I ) 

716 CONTINUE 
721 WRITE (6.608) 

608 FORMAT ( 1H 1 , 22X , 24HTR ANSVE R S E DISPLACEMENTS) 

WRITE (6.609) 

609 FORMAT ( 1H0 , 24X , 20HY INCREASES TO RIGHT/) 

WRITE (6,610) 

610 FORMAT (21X,26HX INCREASES TOWARDS BOTTOM/) 

K= 1 

L=4*N-3 
DO 612 J= 1 , M 

WRITE (6,611) (W(I ) . I =K , L , 4) 

611 FORMAT (8( D11.4,5X) ) 

WRITE (6.613) 

613 FORMAT ( 1H //////////) 

K=L*4 

612 L=L*4*N 

IF (C3.E0.0. ) GO TO 712 
MM=MM*1 

IF (MM.E0.3) GO TO 722 
IF (MM.EQ.4) GO TO 712 

MULTIPLY SINV * 8 * W( NORMAL I ZED ) = WNFXT 

713 DO 705 1=1, INDEX 
SUM2=0. 

DO 704 J=l, INDEX 
$UM1=0. 

00 703 K=1 , I NDEX 

703 SUM1*SUM1+B( J,K)*W(K) 

704 SUM2=SUM2*SINV(I ,J)*SUM1 

705 WNEX T ( I ) = SUM2 



MULTIPLY WNEXT TRANSPOSE * S * WNFXT 



EIGEN 



706 



606 



708 

715 

714 



EIGEN=0. 

DO 706 1=1, INDEX 
DO 706 J*l, INDEX 

E I GEN=E IGEN*WNEXT( I > ♦STIFF ( I , J ) * WNEXT ( J ) 

FIND BUCKLING COEFFICIENT. 

FORMAT ( 1H0,9HITERATI0N,I3,4H, K=,F9.5) 

TEST FOR 20 ITERATIONS. 

IF (L.E0.20) GO TO 710 

TEST FOR CONVERGENCE. 

IF (L.E0.1) SAVE=AK 
IF (L.E0.1) GO TO 708 

IF (0ABS(0ABS(AK)-0ABS(SAVE))/0A9S(AK).LT.. 00001) GO 
SAVE=AK 
L = L* 1 

SAVE W (NORMALIZED) 

DO 715 1=1. INDEX 
WSAVE ( I ) = W ( I ) 

USE RESULTING VECTOR AS NEXT TRIAL VECTOR. 

00 709 1=1, INDEX 



TAKFN FROM 



OF ARRAY) 



IN MAP FORMAT 



TO 711 
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709 W( I )=WNFXT ( I ) 

on t n 7oo 

710 WRITE <*.5 22) 

52 2 Format ( 1H0 , 3 3H20 ITERATIONS WITHOUT CONVFRGf NCF ) 

GC TO 71 2 

711 WRITF (6,575) 

525 FORMAT f 1 HO » 43HB DC K l I NG PARA Mr TF R HAS CONVE RGE 0 TO .001 
MM = 7 

GO TO 722 
717 CONTINUF 
FMO 



ASPtCT RATI 0= 1.00000 POISSON PATIO = 0.0 



NO. ELFMENTS IN X = 4 



NO. 



BOUNDARY CONDITIONS, 1 = FR EF , ?»PINNEI>, 3=CL AMPFD, ^ = SYM. , 5=ANTI-SYM. 
Y=0,fDGF CONDITION 3 Y=B, EDGE CONDITION 4 X=0,F0GE CONDITION 7 > 

RELATIVE LOAD SIZES NX= L.00 N Y= 0.0 NXY= 0.0 



r 8!3V v 6«S%?J3 ST o. l ?9 

-0.09 -0.35 -0.65 0. 72 



-0 . 0° -0.11 0.60 -0.79 

0.22 0.14 0.91 0.14 

0.16 -0.47 0.27 -0.75 

0.39 0.47 -0.98 

0.31 -0.65 -0.57 -0.75 
0.76 -C.33 0.66 -0.68 

0.53 0.34 0.99 0.95 



-0. C9 



O.R 1 ' 

0.71 

0.18 0.58 

-0.55 0.«5 

0.97 -0.60 -0.04 
‘ ““ -0.76 0.97 

-0.19 -0.40 

0.38 — 0 • 3 9 



-O.dl -0. 


74 


0.34 0.99 o.5S -0.31 


I TERATIGN 


1. 


K = 


104.51617 


ITERATION 


7 , 


K = 


26.5343? 


ITERATION 


3 » 


K = 


13.08651 


ITERATION 


4, 


K = 


8.81542 


ITERATION 


5, 


K = 


7.89034 


ITERATION 


6, 


K = 


7. 72830 


I TFRATION 


7, 


Ks 


7. 70094 


ITERATION 


8, 


K- 


7.^9634 


I TFRATICN 


C * 9 


K = 


7.69556 


I TERATICN 


10, 


K = 


7 . 0954,3 


ITERATION 


Ilf 


K = 


7.69541 


31JCKLING PARAMETER HAS CONVERGED TO 





0.68 


-0.54 


0.66 


-0.17 


-0.54 


0.2? 


0.76 


0. 63 


0.19 


-0.17 


0.59 


0. 50 


-0.93 


0.60 


-0.63 


0.17 


0.98 


-0.25 


0.04 


-0.18 


0.30 


0. 15 


-0.61 


0.21 


0.98 


-0.21 


0.21 


-0. 76 


0.55 


-0.33 


-0.97 


0.41 


-0.95 


0.98 


0.75 


0.94 


-0.06 


-0.69 


-0.23 



PERCENT ) 



ELEMENTS IN Y = 4 
= A , FOGE CHMn T T I ON r 
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0.0 



TRANSVERSE DISPLACEMENTS 
Y INCREASES TO RIGHT 
X INCREASES TOWARDS BOTTOM 
0.0 0.0 0.0 



0.0 



0.0 



-0 • 12 12D-01 



-0.3320D-01 



-0 • 50040-01 



-0.56320-01 



0.0 



-0. 1715D-01 



-0.4697D-01 



-0.7081D-01 



-0.79690-01 



0.0 



-0.1214D-01 



-0.3323D-01 



-0.5 009D-01 



-0.5638D-01 



0.0 



0.0 



0.0 



0.0 



0.0 
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